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ABSTRACT 

We study the development of gravitational instability in the strongly non-linear regime. 
For this purpose we use a number of statistical indicators such as filamentary statistics, 
• spectrum of overdense/underdense regions and the void probability function, each of 

Cn| \ which probes a particular aspect of gravitational clustering. We use these statistical 

indicators to discriminate between different approximations to gravitational instability 
which we test against N-body simulations. The approximations which we test are, the 
truncated Zel'dovich approximation (TZ), the adhesion model (AM), and the frozen 
flow (FF) and linear potential (LP) approximations. Of these we find that FF and LP 
break down relatively early, soon after the non-linear length scale exceeds i?* - the 
mean distance between peaks of the gravitational potential. The reason for this break 
' down is easy to understand, particles in FF are constrained to follow the streamlines 

of the initial velocity field. Shell crossing is absent in this case and structure gradually 
. freezes as particles begin to collect near minima of the gravitational potential. In 

LP particles follow the lines of force of the primordial potential, oscillating about its 
minima at late times when the non- linear length scale fe^ L ~ i?*. Unlike FF and LP 
the adhesion model (and to some extent TZ) continues to give accurate results even at 
late times when > i?* . This is because both AM and TZ use the presence of long 
range modes in the gravitational potential to move particles. Thus as long as the initial 
potential has sufficient long range power to initiate large scale coherent motions, TZ 
and AM will remain approximately valid. In relation to AM, TZ suffers from a single 
major drawback - it underestimates the presence of small clumps. Similarly, it predicts 
the right mean density in large voids but misses subcondensations within them. The 
reason for this is clear: The artificial removal of power on scales smaller than fe^ L in the 
■ initial potential in TZ, designed to prevent shell crossing, causes a substantial fraction 

of matter (which would have been clustered in N-body simulations) to lie within low 
density regions at all epochs. On the other hand, TZ is very fast to implement and 
more accurately predicts the location of large objects at late times. 

Key words: Cosmology : theory - galaxies : formation - large scale structure of 
Universe - methods : statistical 
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1 INTRODUCTION 

The Universe on large scales exhibits remarkable struc- 
tural features as demonstrated by the numerous inves- 
tigations of its statistical properties. It is believed that 
this structure arose via amplification, through gravita- 
tional instability, of primordial fluctuations in the den- 
sity of matter. The evolution of such fluctuations can 

© 0000 RAS 



be studied using the well known hydrodynamical equa- 
tions for the gravitating fluid. In the past decade sev- 
eral workers have obtained numerical solutions to these 
equations which confirm that gravitational instability 
can lead to the kind of structure observed in the Uni- 
verse today. 

While numerical N-body simulations are manda- 



tory to approach the precise picture, often our un- 
derstanding of the dynamical processes that lead to 
these structures comes from various approximations 
to the fully nonlinear equations that have been pro- 
pounded. For instance, Zel'dovich showed that gravita- 
tional instability gencrically leads to the formation of 
two-dimensional sheets, the so called pancakes, the ad- 
hesion model which in some sense can be regarded as 
an extension of the Zel'dovich approximation, demon- 
strated that matter moves along pancakes towards fila- 
ments (which form at the intersection of two pancakes) 
and then along filaments towards clumps which form 
at the junction of two filaments (or three pancakes). 
Thus the Zel'dovich approximation and the adhesion 
model showed that gravitational instability leads to 
the formation of cellular structure described by pan- 
cakes, filaments and clumps - a result that has also 
been independently verified by detailed N-body simu- 
lations flZel'dovich 1970| ; fVIelott et al. 198$ |Melott & 
Shand ^rin 1989| ; |Shandarin fc Zel'dovich 1989| ; pVlclott 



Pellma |n fc Shandarin 1994| ; [Bahni fc Coles 1994| ). In re- 
cent years several other approximations to the nonlin- 
ear equations governing gravitational instability have 
been proposed. Such approximation schemes serve a 
dual purpose: Firstly, they have the potential to provide 
us with insight regarding the physical processes which 
led to the formation of structure. Secondly, they are as a 
rule easier to implement and are often computationally 



multiplicity function for overdense regions. It is there- 
fore essential to examine different non-linear approxi- 
mations with a number of distinct (and in some cases 
orthogonal) statistical discriminators. 

In the present paper we compare the following ap- 
proximation methods both with each other and with 
the results of N-body simulations performed using a 
two-dimensional PM code running with 512 2 particles 



on a 512 2 mesh (for details see ( Beacom et al. 1991 )). 
We tested: (a) Zel'dovich approximation (truncated ver- 
sion), (b) the adhesion model, (c) frozen flow approx- 
imation and (d) linear potential approximation. In a 
parallel study, we follow the development of gravita- 
tional instability using a variety of statistical indicators 
which probe its different features. The statistical indi- 
cators which are used for comparing (a) - (d) single 
out certain features of non-linear clustering such as the 
existence of voids, filaments (analogues of pancakes in 
2D) and clumps. The domain of validity of a given ap- 
proximation is therefore discussed with reference to a 
given statistical indicator. 

The treatment followed in this paper extends ear- 
lier work of Coles, Melott & Shandarin (1993) in which 
three nonlinear approximations: the Zel'dovich approx- 
imation, the truncated Zel'dovich approximation, and 
the lognormal approximation were compared with N- 
body simulations. The present treatment is also in 



a sense complementary to recent work (Munshi & 



less ex; )ensive than full N-body simulations. In order to Starobinsky 1994 ; Bernardeau et al. 1994 



Munshi 



apply i, given approximation effectively we should have Sahni & Starobinsky 1994), in which several of the ap- 



a clear understanding of the domain of its validity. It 
might also so happen that certain statistical properties 
are reproduced by an approximation to the same level 
of accuracy as in an N-body simulation although cer- 
tain other statistical properties may be reproduced to 
a much lower accuracy. For instance, a given approxi- 
mation might correctly reproduce the void probability 
function at a given epoch and yet fail to give the correct 



proximations considered by us were examined in the 
weakly nonlinear regime of gravitational instability. A 
common conclusion drawn in the above papers was that 
the Zel'dovich approximation was more accurate than 
either the frozen flow or the linear potential approxi- 
mation when tested against the results of perturbation 
theory in the quasi-linear regime. (This analysis was 
generalised to include Lagrangian perturbation theories 
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in Munshi, Sahni & Starobinsky (1994). The present pa- 
per presents a fully nonlinear treatment of the problem 
thereby considerably extending the quasi-linear analy- 
sis of the above authors. Our tests were conducted in 
two dimensions and complement the three dimensional 
analysis of Melott et al. (1994) and Melott, Shandarin 
& Weinberg (1994). Some of our results may however 
be carried over to three dimensions as well. The exact 
formulation of the truncated Zel'dovich approximation 
used here may be found in Melott, Pellman & Shandarin 
(1994). 

The paper is organised as follows. In section II we 
briefly discuss the various approximations to gravita- 
tional instability that we have chosen to compare with 
N-body simulations. Section III is divided into several 
subsections each of which deals with a particular statis- 
tic. In section IV we summarize the chief results of our 
investigations. 

2 NONLINEAR APPROXIMATIONS 

Consider pressureless matter with density p(t,x). The 
dynamics of such a fluid is governed by the expansion of 
the Universe as also by inhomogeneities in its distribu- 
tion. The component of the velocity which arises solely 
due to inhomogeneities in the density field is known 
as the peculiar velocity v(t,x) = r — Hr. The com- 
bined evolution of the density p, peculiar velocity and 
the peculiar gravitational potential (p(t, x) is given by 
the following well known system of coupled nonlinear 
equations: 



^ + 3Hp+-V-(pv) = 0, 
at a 

dv 1 1 
— + -(v V)v + tfv = — Vip, 
at a a 

\7 2 ip = AirGa 2 (p-p ), 



(1) 

(2) 
(3) 



where a(t) is the cosmic expansion factor, H is the Hub- 
ble parameter and po is the average density of the fluid. 
The spatial derivatives in (||) - (0) are defined with re- 



— + u-V u = -— {u + A\/ip ), 
aa 2a 



spect to the comoving coordinate x = r/a. Choosing a 
new time variable a(t), and defining a comoving velocity 
variable 

G^X V /A\ 

u = -r = — > 4 

da aa 

we obtain the following form for the Euler equation (^|) : 

(5) 

where A = 2/(3H 2 a 3 ) is a constant for a flat Universe 
with dust-like matter. In cosmological problems, where 
initial perturbations correspond to the growing scalar 
mode, the velocity field u is potential u = — V$ until 
multistream regions develop. 

At earlier moments of time when inhomogeneities 
are small the solutions to equations ([]])- (||) may be ob- 
tained by linearization. We then have during the linear 
stage 



u(q) = -AV<p(q) 



(6) 



where q are the initial (Lagrangian) coordinates. There- 
fore, initally, the velocity potential and the gravita- 
tional potential are simply proportional to one-another 
$(q) = Aip(q), both unchanging in a flat matter domi- 
nated Universe. 

Later, in the nonlinear regime, there is no easy solu- 
tion to the basic system (Q)-@. The different nonlinear 
approximations considered by us can be conveniently 
described as different ways of simplifying equation (||). 

2.1 Truncated Zel'dovich Approximation (TZ) 

The Zel'dovich approximation (henceforth ZA) may be 
obtained from (||) by setting its right hand side to zero: 



Du du , _. 

— = — + u-V u = 0. 

Da aa 



(7) 



where D/Da is the convective derivative. The above 
equation says that the dynamics of the fluid element 
is governed purely by "inertia" . It has an immediate 
solution in terms of the displacement of the fluid el- 
ement from its initial position with constant velocity 



dZel'dovich 197C| ) 
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x = q + a(£)u(q) 



(8) 



By setting the right-hand side of eq.(Q) to zero 
in ZA, one extrapolates the linear relation (^|) be- 
tween velocity and gravitational potential into nonlin- 
ear regime where the potentials are generally time de- 



pender lt u(x. t) — — jlVy(x, t). 



ZA works reasonably well so long as streamlines of 
flows do not cross one another. However, multistream 
flows invariably form at the locations of pancakes, which 
grow progressively thicker leading to the ultimate break 



down of the Zel'dovich approximation (Shandarin & 
Zel'do vich 1989| ). An extension of ZA called the trun- 



cated Zel'dovich approximation (Coles, Melott & Shan- 



darin ; 993) is based on the observation that the for- 
mation and thickening of pancakes can be delayed by 
artificially removing power on all scales smaller than 
the one that is currently going nonlinear. The length 
scale and window shape with which the original spec- 
trum is best smoothed has been determined in three 
dimensions by Melott, Penman & Shandarin (1994). In 
our simulations we have used a fc-space Gaussian win- 
dow exp(— fc 2 /2fcg) to implement the necessary trunca- 
tion. As found by Melott, Penman & Shandarin (1994). 
the optimal cutoff scale is related to the scale entering 
nonlinearity. However, the precise filtering scale fc^ 1 de- 
pends on the spectrum. This represents a drawback of 
the model as we cannot obtain the best results for an 
arbitrary spectrum based on first principles. However, 
the cutoff scale is only weakly spectrum dependent. The 
optimal value of fcg 1 for spectra considered in this pa- 
per will be discussed in the next section (cf. Table II). 
Nevertherless, let us note that even if we do not use the 
best filter for a given spectrum but a fixed spectrum in- 
dependent filtering k(j ^ kcin) the approximation still 
retains most of its positive features. One major advan- 
tage of the TZ is the extreme simplicity of its imple- 
mentation. 



An extension of TZ is the use of second order per- 
turbation theory combined with the smoothing of the 
initial potential. This produces somewhat more accu- 
rate but not qualitatively different results from TZ 



(Buchert, Melott k Weiss 1994; Melott, Buchert, and 



Weiss 1994). We do not include this second-order ap- 



proach in this study. 

2.2 Adhesion Model (AM) 

The adhesion model is an extension of the Zel'dovich 
approximation. In the adhesion approximation the right 
hand side of (||) is replaced by an artificial viscosity 
term to mimic the effects of nonlinear gravity on small 
scales and to stabilise the thickness of pancakes. The 
resulting equation is the well known Burger's equation 



and has the form (Burgers 1974; Gurbatov, Saichev, & 



Shandarin 1985,1989) 



du 

da 



+ (u.V)u = i/V 2 u, 



(9) 



where v is the coefficient of viscosity. It is interesting 
that in the limit v — > the right hand side of (^) 
remains finite only in those regions where large gradi- 
ents in the velocity field exist (viz inside the pancakes) , 
and vanishes elsewhere. As a result the adhesion model 
reproduces the results of the Zel'dovich approximation 
exactly in regions outside of the pancakes themselves. 
Accordingly, the adhesion model reduces to ZA for the 
early time moments or sufficiently smoothed initial con- 
ditions when no shell-crossing is present. For vanishing 
v the adhesion model has an elegant geometrical inter- 
pretation which we have used to construct the skele- 
ton of the large scale structure predicted by this model 
QGurbatov, Saichev, fc Shandarin 1985,1989| ; |Po gosyan| 



1989; Kofman, Pogosyan, & Shandarin 1990; Kofman 



et al. 1992; 3ahni, Sathyaprakash & Shandarin 1994) 



An undesirable limitation of the geometrical prescrip- 
tion is that it does not give particle positions but only 
locations of filaments and clumps which have to be 
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smoothed by an appropriate filter in order to lend them- 
selves to a comparative treatment with other models 
and with N-body simulations. A study using particles 



dynamics of the linear potential (or frozen potential) 



approximation (Brainerd, Scherrer & Villumsen 1993 



^=0 
oa 



(Melott, Shandarin, & Weinberg 1994) in three dimen- 
sions shows general agreement with our results when 
equivalent tests were done. 

2.3 Frozen Flow Approximation (FF) 

The underlying philosophy of FF is in a sense just the 
converse of ZA since the inertia of particles is neglected 
in this approximation which requires particles to con- 
stantly upgrade their velocity to a value determined by 
the local value of the linear velocity field. More pre- 
cisely, FF corresponds to neglecting both the nonlinear 
term, namely, (u • V)u and the right-hand side in (g) 
( [Matarrese et al. 1992; ) 

(10) 

so that the comoving velocity field remains fixed to its 
linear value u(x, t) — u(q = x). It is clear that matter 
in FF is collected with time in the points u(q) = i.e. 
in the positions of the local minima of the initial gravi- 
tational potential. Therefore FF cannot be expected to 
work even qualitatively for late time moments when the 
scale of nonlinearity escalates above the typical distance 
between minima of the initial potential. 

2.4 Linear Potential Approximation (LP) 

N-body simulations show that the gravitational poten- 
tial evolves much more slowly than the density field 
( [Brainerd, Scherrer fc Villumsen 1993| ). This is so be- 
cause relative to 5 the potential tp is dominated by 
small-fc modes which obey the precepts of linear theory 
longer than large-fc modes. The linearized Poisson equa- 
tion V 2 ^ = AnGa 2 8p demonstrates that p ~ const, as 
long as S -C 1 and the Universe is flat and matter domi- 
nated. Extending this assumption (tp ~ const.) into the 
nonlinear regime as well, we arrive at the following gen- 
eralisation of the Euler equation which describes the 



Bagla & Padmanabhan 1994) 



9u , _, 3 , . ¥ . 

— + u • V u = - — (u + AV<p ), 
oa 2a 



(11) 



where tp = tp(x., to) — <£>(q). This equation defines the 
force acting on a fluid element at the instant a(t) using 
the primordial value of the potential (/> = i^o- In a sense 
the LP can be regarded as an N-body simulation in 
which the value of the potential is not upgraded after 
each time step. 

Both TZ and AM have single-step analytical solu- 
tions. Consequently, there is no need to evolve the fluid 
iteratively; given some initial conditions these approxi- 
mations have the ability to directly give the configura- 
tion at any epoch which may be of interest. In contrast, 
LP and FF are Eulerian approximations and have no 
analytical solutions except in some special cases. Op- 
erationally they are similar to full N-body simulations 
which evolve the fluid iteratively, except for the fact 
that in LP the potential is kept frozen to its initial value 
and in FF neither is the velocity potential upgraded nor 
is particle inertia taken into account. Since PM type N- 
body simulations are easy to do on modern computers, 
it is not clear whether LP and FF have value beyond 
the descriptive insight which they provide. 

3 COMPARATIVE STUDY OF VARIOUS 
APPROXIMATIONS 

In this section we employ a number of statistical tools 
to compare the approximations mentioned in the previ- 
ous section with N-body simulations. We use the same 
initial conditions for all the approximations and they 
form a subset of the initial conditions used by Beacom 
et al. (1991) and Kofman et al. (1992) for other pur- 
poses. Time evolution of the N-body models can be 
seen in the video accompanying Kofman et al. (1992). 
All our comparisons are carried out in two dimensions 
with the initial potential being specified on a grid of 
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size 512 x 512. More specifically, the models for which 
we have carried out the comparison are either feature- 
less or truncated power law spectra of the general form 



, , k n , for k < k c 
^ W K l 0, for k > k c . 



(12) 



We have considered three different spectral indices, n — 
2, 0, — 2, with a cutoff in each case, at the Nyquist 
wavenumber: k c = 256k t, where kf is the fundamental 
mode. In addition to this, we have ati = model with 
a truncation at k c — 32k f which serves to illustrate 
the effect of an abrupt cutoff in the power spectrum as 
happens in some models of dark matter like hot dark 
matter. Thus, we have a total of four models in all. 

All our simulations of the various approximations 
are performed using a particle code excepting the ad- 
hesion model which is simulated using the well known 
geometrical interpretation of the solution to Burger's 



equation (Shandarin & Zel'dovich 1989; Pogosyan 



1989; ;ahni, Sathyaprakash & Shandarin 1994). Conse 



quently, we could not include adhesion in studying some 
statistical properties, such as the position correlation 
coefficient or filament statistics, which rely on having 
particle positions. Where we could compare AM with N- 
body we expect the former to do somewhat better than 



what our results convey. See also ( Melott, Shandarin, fc 



Weinb erg 1994 ) for a particle-based three-dimensional 
study of AM. 

We compare the evolved density fields and the 
quantities derived from them when different scales are 
going nonlinear. We choose ct(/cnl), the epoch when 
the scale 2ir/k^ is going nonlinear, as a convenient 
measure of "time" with which to characterize different 



(13) 



regimes in nonlinear gravitational clustering: 

( ; fc fc ; p(k)kdk \ 1/2 

o-(fcNL) = — r 

{j^P(k)kdk) 

Here &n is either the Nyquist wavenumber or the cutoff 
mode k c , whichever is smaller. For truncated power law 
spectra (with a cutoff at k c ) 



cr(k NL ) = - — 

V fc NL 

and 



. In (fe c /A;/) 
owl) = I . I ■ " 



-2. 



(14) 



(15) 



ln^NL/A;/) 

The first scale to go nonlinear is the one corresponding 
to either the Nyquist wavenumber (in the case of mod- 
els with no cutoff) or the mode k c . When this happens, 
by definition, a — 1. As a increases, successively larger 
scales enter the nonlinear regime. For concreteness we 
have considered in our comparison those values of a for 
which the scales going nonlinear are k c , k c /2, etc., and 
we stop the integration when the scale entering the non- 
linear regime comes close to the size of the simulation 
box. 

In this study we suggest that two natural scales 
characterizing a given model may be well suited for giv- 
ing bounds on the validity of some approximation meth- 
ods. These are: (i) the scale i?*, corresponding to the 
average distance between the peaks of the potential and 
(ii) the scale R v characterizing the correlation length of 
the potential. They are given in terms of the moments 
of the potential field by the following expressions: 



R = V2^ R* = V2^ 
where the moments <jj are defined by 



(16) 



a] cx / k 2j ~ i P{k)kdk. 



(17) 

The epoch cr* (cr v ) when the scale R* {Rip) is going non- 
linear can be found from (|lj) and ( |l5| ) by substituting 
&nl = R* 1 {R^ 1 )- The values of er* and cr^ are listed 
in Table I for the various models under discussion. The 
values of i?* and R v as well as 7 = R*/R v are plotted 
in Fig. [l] as functions of the spectral index n. 

The values of i?* and R v as well as 7 = are 
plotted in Fig. [l] as functions of the spectral index n. 

We recognize that the specific values of ^(cr^) and 
R* {Rip) are often (and in particular in our case of power- 
law initial potentials) determined mainly by numerical 
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k c = 32k f k c = 256fc^ 



11 


(7* 




£7* 


(Tip 


-2 


2.36 


oo 


2.55 


CO 





3.75 


20.3 


4.71 


120 


+2 


4.00 


14.0 


4.00 


22.2 



Table 1. The scales R* and R v of potential corresponding to the average distance between the peaks of the potential and the correlation 
length, respectively. The box is assumed to be of unit length. Also tabulated are the corresponding epochs ex* and a v when these scales 
go nonlinear. 



cutoffs introduced by the limitations of our computer 
simulation. In fact this reflects once again the effect of 
the finite grid used in any simulations on the represen- 
tation of the underlying initial spectrum. For the real 
Universe and spectra such as CDM, physical cutoffs are 
provided by the horizon scale (for R v ) and by the free- 
streaming distance (for i?*). 

As mentioned earlier in Sec. [2.l| the optimal 
smoothing scale that needs to be used in TZ simulations 
depends on the spectrum. By definition, the optimal 
smoothing scale is that scale which obtains the max- 
imum correlation coefficient of TZ density fields with 
N-body density fields. We have found that its relation 
to the scale entering nonlinearity is fairly independent 
of the epoch. (If this were not so then the very con- 
cept of truncated Zel'dovich approximation would lose 
its meaning.) Table II lists the optimal smoothing scales 
fcg 1 = k~ pi for different spectra considered by us. 

3.1 Visual comparison 

To begin with, we make a visual comparison of the var- 
ious approximation schemes with N-body simulations. 
The structure obtained using the AM and the evolved 
particle positions in the case of FF, LP and TZ are 
shown in Fig. ^|a-d for four different spectra. From top 
to bottom, the pictures correspond to N-body, adhesion 
model, frozen-flow approximation, linear potential ap- 
proximation, and truncated Zel'dovich approximation, 
respectively. Fig. |^a corresponds to n = 0, k c = 32k / 
model and Fig. &> to n — 0, k c — 256k f model. In Fig. 



||c and ||d the spectral index is n = 2, and —2, respec- 
tively, and k c — 256/c / in both the cases. In Fig. ||a, b 
and d the left hand panels correspond to an epoch a 
such that a ~ <7* (k^ ~ i?*) and the right panels to an 
epoch a ~ a v (k^l ~ R^). In Fig. ||c both the left and 
the right hand panels correspond to an epoch a > a v . 
(In the case of n = 2 models there is a lot of small scale 
power. Consequently, the pictures at a stage when i?„ 
is going nonlinear looks too grainy. Therefore it is not 
easy to compare them visually at that epoch.) 

We find that at the epoch when the scale going 
nonlinear is R„, and at earlier epochs, all the approxi- 
mation schemes appear to reproduce the structure with 
roughly the same accuracy as in N-body simulations. 
(This is also reflected by the high value of the corre- 
lation coefficient before the epoch <r* as discussed in 
sections |3.2.l| and |3.2.2| - see Fig. | and g). 

The epoch corresponding to fc^ L ~ i?* characterizes 
the completion of cellular structure which forms from 
an initially smooth distribution of matter. Later epochs 
are characterised by the relative motion and mergers of 
the structure elements governed by mutual attraction of 
large mass concentrations (knots and filaments) as well 
as repulsion from underdense interiors of cells (voids). It 
is clear that both FF and LP which fix the structure on 
scale i?* are unable to describe this process even qual- 
itatively, and therefore begin to fail beyond the epoch 
it*. These conclusions are borne out by Fig. |2|a-d. 

From the right panels it is clear that close to the 
epoch when the scale going nonlinear is R v the struc- 
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11 


fcopt/fcNL 


-2 


0.5—1.50 





1.25 


+2 


1.00 



Table 2. The optimal smoothing scale fc * used in the truncated Zel'dovich approximation scheme depends on the index of the power 
spectrum as shown above. For n = —2 spectra the density correlation coefficient is virtually the same for a wide range of values of fc op t- 



ture ot tained using AM is still in excellent visual agree- 
ment with that of N-body simulations. TZ has a reason- 
ably good visual agreement with N-body on large scales 
though at this epoch the small scale features abun- 
dant in N-body simulations (especially for spectra with 
n > 0) are not present in the TZ simulation. However, 
much before this epoch, in fact soon after the scale R„ 
has crossed nonlinearity, FF and LP approximations fail 
to give the right picture. Particles in FF, approach the 
valleys and the minima of the potential asymptotically, 
leading to greatly thinned out filaments which eventu- 
ally empty out and vanish as the particles gradually fall 
into the minima of the potential. In the case of LP, par- 
ticles execute oscillations around the troughs of the pri- 
mordial potential partially simulating, at earlier epochs, 
the results of N-body simulations wherein the pancakes 
neither thicken as in the Zel'dovich approximation nor 
do they thin out as in FF. The contour diagrams of the 
initial potentials used in our simulations, shown in Fig. 
|3|, bear out this claim about the behaviour of particles in 
FF and LP. Notice especially the pictures correspond- 
ing to the n — 2 power-law model (cf. Fig. ^|c) which 
has a lot of small scale power. Here we see that the 
particle positions in FF and LP are essentially frozen 
beyond the epoch fcj^ L ~ i?*. While the dynamics in 
LP and FF at all times are determined by the gradients 
of the local primordial potential, we know from N-body 
simulations that beyond the epoch of formation of cellu- 
lar structure, the small scale features of the primordial 
potential play little, if any, role in the furtherance of 



gravitational clustering ( Beacom et al. 1991 ; Pogosyan 
199C |; Rittle, Weinberg, Park 199~| |Pauls and Melott 



1994). It is therefore to be expected that LP and FF 
will not be able to reproduce qualitatively the features 
of hierarchical clustering. 

The visual agreement of the pictures obtained using 
AM and TZ with N-body lasts for epochs much longer 
than that for cither FF or LP since the former two ap- 
proximations successively use power on larger scales to 
influence the dynamics of the fluid. Even though the 
local gravitational potential has changed substantially 
by the time i?* has gone nonlinear, it has not evolved 
so much as to compete with the effect of power on large 
scales. Thus, as long as the initial potential has suffi- 
cient large scale power to give rise to coherent motion 
over large scales, TZ and AM will remain approximately 
valid. Following Kofman et al. (1992) we speculate that 
the coherence length R v of the primordial potential is 
important in this discussion. By the time a v when the 
scale Rtp goes nonlinear, in fact even at slightly earlier 
epochs, both AM and TZ begin to produce structure 
substantially different in small detail from that seen in 
N-body simulations. However, it is impossible to deter- 
mine from our results whether this change is due to a 
transition at R v or our increasing ability to resolve de- 
tail as the simulation goes nonlinear on larger scales. 
TZ does not produce small objects, and AM puts them 
in the wrong place. Pauls & Melott (1994) have shown 
that even at much later times than a v , the primordial 
potential can determine the coherent motion of large 
chimps, and TZ can produce correct positions for them 
while AM begins to make errors in the position of large 
objects as well and all other approximations have bro- 
ken down long before. However, let us stress that both 



8 



TZ and AM continue to reproduce qualitative features 
of the structure even beyond R v . 

The TZ approximation suffers from one major 
drawback. Since in this approximation we have arti- 
ficially removed power on small scales, matter never 
gets collected in small clumps. It does, however, put 
about the right amount of mass in large clumps. Conse- 
quently, a substantial fraction of matter (which would 
be in small clumps in an N-body simulation) lies within 
low density regions at all epochs. As a result TZ (which 
has the advantage of being computationally very fast) 
is not well suited for studying small clumps even though 
it gives a remarkably good correlation coefficient when 
compared with N-body simulations. Similarly TZ gives 
the right mean density in large voids but misses subcon- 
densations within them. This is the price it pays for its 
greater accuracy in locating the large-scale mass distri- 
bution. These views, based largely on a visual compar- 
ison of the different approximations, are borne out by 
more quantitative comparisons which we discuss below. 
The AM suffers from two major drawbacks. First, it is 
computationally expensive, sometimes approaching the 
cost of an N-body simulation. Second, although it does 
produce small clumps, it occasionally makes major er- 
rors in their position (sometimes comparable with the 
scale of nonlinearity), especially at late times and larger 
n. 

3.2 Quantitative comparison 

We now turn to a quantitative comparison of the various 
approximation schemes with N-body simulations. The 
information about any statistic can in principle be in- 
ferred from a knowledge of all the N-point correlations, 
or at least a substantial number of them. In practice, 
however, the full hierarchy of correlation functions is 
virtually impossible to calculate due to heavy comput- 
ing requirements. Even if it were in principle possible to 
compute the higher order correlation functions it is not 



clear how much light that would shed on structural units 
that are often of interest such as clumps and filaments. 
It is with the aim of studying such structural units that 
we have chosen a number of statistical tools each of 
which addresses a specific structural feature present in 
the simulation. To this end we choose the time evolution 
of the following indicators for comparison: 

(i) Correlation coefficient of particle positions of dif- 
ferent approximation schemes with N-body, 

(ii) Correlation coefficient of density fields of different 
approximation schemes with N-body, 

(iii) Number of clumps, (defined as regions of a cer- 
tain overdensity) , 

(iv) Filament statistics, 

(v) Number of voids, (defined as regions of a certain 
undcrdensity), 

and 

(vi) Void probability function. 

We note that the first two indicators are primarily 
dynamical, in that they test for specific point-by-point 
agreement between the approximation and the N-body 
simulation. The others are global statistics, and test for 
particular kinds of similarity. 

3.2.1 Correlation coefficient of particle positions 

Given an approximation scheme an obvious question 
that comes to mind is how well can the approximation 
scheme evolve the particles in relation to the exact N- 
body simulations. To answer this question we consider 
the correlation coefficient of the coordinate positions 
of an approximation scheme and N-body. Let X\(er) 
and X^(cr) denote the position of the i th particle, at 
an epoch a(t), in an approximation A and N-body, re- 
spectively. The displacement of the particles from their 
unperturbed coordinates X(cr(fo)) is given by: 

AX( ff (t)) = X(a(t)) - X(a(t )) (18) 
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The linear correlation coefficient of two vector fields 
AX A (<j(i)) and AX N (cr(i)) is defined by 

x , — zm^m (19) 

where a dot denotes the scalar product of the two vector 
fields, where <5X = AX — (AX) is the deviation of the 
displacement vector from average displacement, where a 
summation is over the entire sample and where ( ) indi- 
cates average over the entire sample. This is a straight- 
forward generalization of the familiar correlation coeffi- 
cient defined for scalar fields. In Fig. |i| we have plotted 
the evolution of rx for different approximation schemes, 
and for different spectra, as a function of cr(t). Here, 
and in Fig. |] and [To], we adopt the following scheme 
for displaying the evolution of different statistics: The 
top left panel corresponds to n = 0, k c = 32k f model 
and the top right panel to n — 0, k c = 256k f model. 
The bottom left panel corresponds to n — 2 model and 
the bottom right panel to n = —2. In Fig. [| the re- 
sults AM are missing since our implementation of this 
approximation scheme (using an osculating paraboloid) 
does not obtain particle positions. 

We observe that when the spectral index n = — 2 
there is an excellent agreement between all the three 
approximations and the N-body. (rAx is always larger 
than about 0.9 for all approximations for this spec- 
trum.) In the n = 2 case, surprisingly, the correlation 
vanishes to begin with, building up later to reach a max- 
imum of about 0.6 before dropping. We do not under- 
stand this behaviour entirely but suspect that some nu- 
merical effect from the excessively large power on very 
small scales present in this case is the root cause for 
such a behaviour. In this case we see that none of the 
approximations produce the right kind of displacement 
of particles. This is absolutely in agreement with Fig. 

(corresponding to n = 2 spectrum) wherein we see 
that in the case of both FF and LP matter simply gets 
collected into the local minima of the potential without 



ever transferring the power to larger scales while in the 
case of TZ matter does not cluster on small scales at 
all. Although the agreement of all approximations with 
N-body is relatively poor for this spectrum, we find that 
TZ gives consistently higher values for the correlation 
coefficient especially at late times when clustering is 
more prominent on large scales. 

From Fig. |J we find that for spectra with substantial 
power on a wide range of scale, such as n = or n = 
2, particle positions in TZ agree with those in N-body 
much better than do either LP or FF. This is because 
of the fact that for such spectra, wherein both small 
and large scales dictate the dynamics (with larger scales 
being more important at later epochs) LP and FF break 
down at relatively earlier epochs than TZ and AM. 

3.2.2 Correlation coefficient of density fields 

While the correlation coefficient of particle positions 
tells us precisely how the structure is produced in an 
approximation scheme it is seldom the main quantity 
of interest; it is only a measure of how good a dynam- 
ical approximation scheme is in relation to the N-body 
solutions. One is often interested in the evolution of 
the density field since it lets us infer the evolution of 
many other structural units such as clumps, filaments 
and voids. In order to study the time evolution of the 
density field we obtain the density field by employing 
the cloud-in-cell (CIC) algorithm. This algorithm can 
only be used when particle positions are known and 
hence cannot be directly used for our AM simulations. 
In the latter case, we use the structural units of clumps 
and filaments (plus the "free" particles - those that have 
not yet fallen into caustics ) given by the model to re- 
construct the density field. The density field for AM 
simulations is obtained in three steps: (i) The mass in 
the "free" particles is distributed using the CIC algo- 
rithm, (ii) Each clump given by AM is assumed to be a 
"Gaussian hill" with the variance of the Gaussian cho- 
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sen to be proportional to its mass, (iii) The rest of the 
mass is distributed uniformly amongst filaments which 
is then smeared using the CIC algorithm. 

The density fields so obtained are smoothed at a 
certain scale before computing the correlation coeffi- 
cient. Such a smoothing is motivated by the fact that 
the density fields evolved by approximation schemes arc 
not expected to agree with those of N-body simulations 



and LP. For the n = 2 model AM, FF and LP, all give 
roughly the same correlation at all epochs, while TZ is 
better. In contrast to the pictures wherein the agree- 
ment of TZ with N-body is not so remarkable, the den- 
sity correlation is extremely good up to very late times. 
In absolute terms however, both TZ and AM eventually 
break down. The results of AM shown here are entirely 



consistent with those found before ( Mclott, Shandarin 



in greet detail; any agreement is to be expected only fc Weinberg 1994) using a particle AM code, in spite 



after the small scale inhomogeneities are smoothed out. 
In fact, the relevant question here is: "How well can a 
given approximation mimic the results of exact equa- 
tions on medium scales?" . The physical reason for this 
is that the very large scale properties of the universe can 
be studied by ZA or Eulerian perturbation theory, while 
the smallest ones by N-body simulations plus hydrody- 
namics. We might mention that an improvement in the 
correlation coefficient for TZ is expected if instead of 
filtering the density on a fixed scale, say, k ~ 64fcy, a 
variable filter scale which is a constant multiple of the 
scale of non-linearity k oc /cnl is chosen. This is demon- 
strated in Fig j^b for a density which is filtered on a 
scale k = 2/cnl for the spectrum n = 0, k c = 256k f. 

Following Coles Melott & Shandarin (1993) we use 
the correlation coefficient of density fields to compare 
the approximation schemes with N-body simulations. 
Given density fields pa(c) and pn(c), corresponding to 
an approximation scheme A and the N-body simulation, 
respectively, the correlation coefficient of these fields is 
defined by a formula similar to equation (p"9|): 



rs 



£i(*i) 3 E*(*£) 2 



1/2 



(20) 



where 8 = (p — Po)/Po is the density contrast. 

The evolution of the statistic rg is shown in Fig. 
|B|a. The arrangement of the panels here is as in Fig. |I| 
Here we have also included the results of the adhesion 
model. We notice that except for the n = 2 model TZ 
and AM are in better agreement with N-body than FF 



of the artificial smoothing of filaments and clumps that 
has gone into our code to make the geometrical method 
mimic the particle method. In the case of n = power- 
law models, we note that soon after the epoch <t», say 
cr = 8, while the LP has managed to produce a correla- 
tion coefficient of about 0.5, the FF has a far less value 
at this epoch. However, both the AM and the TZ give 
values substantially larger than 0.5 indicating their va- 
lidity at this and later epochs. The reason why the FF 
does so badly can be traced to Fig. ||a-d wherein we see 
that matter has completely emptied out into rivulets of 
the potential wells by the epoch er* . 

The superior performance of TZ on scales down to 
^nl combined with its speed makes it a good approx- 
imation for studying the mass distribution from scales 
of about £> + 10 13 Mg) on up, where b is the bias 
parameter and z is the redshift. As we shall see, it fails 
on smaller scales. But the mass resolution is two orders 
of magnitude better than in the past. 

A good density correlation of an approximation 
scheme with an N-body simulation does not necessarily 
mean that other statistical indicators too will give good 
results when testing the approximation with N-body 
(unless of course the density correlation approaches 
unity). Conversely, a poor density correlation does not 
necessarily imply that the approximation will also pre- 
dict incorrect results for other statistics. However, our 
two previous statistics, the position and density corre- 
lations do test the dynamical accuracy of our approxi- 
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mations. The tests which follow probe agreement with 
respect to global statistical measures, which is a some- 
what different issue. Correct global statistics do not nec- 
essarily imply correct dynamics. In what follows, there- 
fore, we supplement the two correlation coefficients dis- 
cussed above with other statistical indicators that ad- 
dress issues relating to the ensemble of structural units 
present in the simulation. They do not test whether 
specific individual units are in same the same place but 
they do test an overall "resemblance." 

3.2.3 Evolution of the number of clumps 

An important statistic which any theory of LSS is ex- 
pected to explain is the mass function of galaxies or of 
clusters of galaxies: how the total mass is distributed in 
objects of different masses. A related question is how the 
number of clumps, defined as connected regions of a cer- 
tain overdensity, evolves with time. Due to lack of space, 
here we address only the latter, relegating the more im- 
portant former question to a future work. The num- 
ber of clumps, at any moment, clearly depends on the 
density threshold p c we use to identify regions of over- 
density. However, since the aim here is to compare the 
predictions of different approximations with N-body, it 



evolution of the number of clumps and the bottom pan- 
els show the evolution of the fraction of mass in clumps. 
In Fig. |?]a left panels correspond to n — 0, k c — 2,2k / 
model and the right panels to n = 0, k c — 256fc/ model. 
In Fig. 0b the left panels correspond to n — 2, and the 
right panels ton = —2, models, respectively. The results 
of N-body are shown in thick solid lines. Following the 
N-body curve we see that generically, there are two dis- 
tinct phases in the clustering of matter via gravitational 
instability. During the first phase the number of clumps 
keeps increasing reaching a maximum after the epoch 
er„ at which time the formation of cellular structure is 
complete. By this epoch nearly 50 % of the matter that 
ever gets bound has fallen into the local wells of the ini- 
tial potential causing major changes in the local value 
of the potential without, however, disrupting large scale 
modes. During the second phase clustering proceeds hi- 
erarchically, with smaller clumps merging with one an- 
other to form clumps of larger mass. As a result the 
density contrast of the clumps alone (not shown) keeps 
building up at a phenomenal rate whereas the number 
of clumps begins to fall. Roughly, c* characterizes the 
epoch of the transition from the cellular to the hierarchi- 



cal phase of clustering ( Sahni, Sathyaprakash fc Shan- 



hardly [matters what density threshold we choose pro- darin 1994|) . The gravitational potential on small scales 



vided it suffices to obtain well defined clumps. In Fig. 
^ we have shown regions in our N-body simulations of 
density p > p c at two different epochs each, for the 
models n = 0, k c — 32fc/, and n = 0, k c = 256/c/. 
We see that the clumps are well defined for the cho- 
sen density threshold and we use appropriate density 
thresholds for the different spectra considered by us. 
A clump is now defined as a connected region, in the 
sense of a "friends-of-friends" algorithm, of overdensity 
greater than or equal to the density threshold p c . 

In Fig. ^ we have shown the evolution of the number 
of clumps for different simulations of the various power- 
law models discussed earlier. The top panels show the 



undergoes substantial changes during the second phase 
culminating in the disruption of any initial small-scale 
coherence that might have existed in the primordial po- 
tential. Further discussion of the evolution of the poten- 
tial can be found in Pauls & Melott (1994). Beyond the 
epoch cr* no approximation scheme which proposes to 
keep the local values of the potential unchanged, and 
at the same time does not use large scale power in de- 
scribing the dynamics, can predict correct gravitational 
clustering. 

Strictly speaking, as far the evolution of clumps is 
concerned the agreement between N-body simulations 
and approximation schemes, depends on the spectral 



12 



index. However, based on Fig. [7] we first make the fol- 
lowing general remarks. We note that with the excep- 
tion of AM none of the other approximations reproduce 
the expected fall in the number of clumps at late times 
caused by merger of clumps. In FF and LP the num- 
ber of clumps, at earlier epochs, show the general trend 
of a sharp rise and closely resemble the predictions of 
N-body simulations till about cr*. However, neither of 
these approximations show a proper fall off in the num- 
ber of clumps. In fact, with the exception of the n = 0, 
k c = 32k f model, there is no fall off seen in the num- 
ber of clumps in these approximations which is a major 
feature of gravitational clustering: In these approxima- 
tions there is relatively little evolution in the number 
of clumps at later times. This is in contrast to N-body 
simulations wherein the merger of clumps, with smaller 
clumps falling into the wells created by the larger ones, 
is a never-ending process. The bottom panels of Fig. 

lend further support to the viewpoint that these two 
approximations do not predict the correct clustering of 
matter beyond the epoch <r*. We observe that in N- 
body simulations matter gets continuously drained into 
clumps whereas in none of the approximation schemes, 
except in AM, is this phenomenon seen generically. We 
note that the statistics of clumps produced by TZ is 
almost never in agreement with N-body simulations. 
The reason for this is that in this approximation the 
clumps are at no time well defined objects. Clumps 
are relatively short scale features which form soon after 
the rms linear density contrast reaches unity. There- 
after they mature, and gain identity by the epoch when 
the scale going nonlinear is However, in the case 
of TZ the linear theory rms density contrast is never 
allowed to greatly exceed unity: power on successively 
larger scales, in fact power on roughly the scale that 
is entering nonlinearity, is filtered out. Consequently, 
although voids are well defined in this approximation, 
clumps never acquire a permanent identity. The bot- 



tom panels of Fig. [7|a show that hardly 50 % of the 
total mass ever gets collected in clumps in TZ. This, in 
addition to the very low number of clumps that the ap- 
proximation predicts, makes the approximation scheme 
unsuitable for any study concerning small scale features 
of LSS such as individual galaxies. This is a shortcom- 
ing of this approximation, but one that is not entirely 
unexpected. 

In our opinion the AM is best suited in describing 
clustering which statistically resembles N-body simu- 
lations. It predicts the right kind of growth law for 
the fraction of matter in clumps as well as the right 
merger histories of collapsed objects (except perhaps in 
the n = — 2 case). 

3.2.4 Filamentary statistics 

The second in our study of the structural units of LSS 
is filaments. In order to understand the formation and 
evolution of filaments we employ a statistic first sug- 
gested by Vishniac (1979) and later used by Nusser and 
Dekel (1990) to study filamentarity in different models 
of structure formation. This statistic is obtained by first 
identifying the moments of the distribution of particles 
around a chosen centre and then constructing a scalar 
from these moments. 

Let Mg(R) and (i?) denote the first and sec- 
ond "moments" of the distribution of the particles lo- 
cated within a distance R around the k th particle: 

-, N k 

M *=^E^> (21) 
<" = ^I>;4 ( 22 ) 

Here is the number of particles within a distance 
R from a centre located at the k th particle, and xf, 
a = 1,2, denotes the position vector of the j th particle 
relative to the k th particle. The filamentary statistic 
S(R) is the ensemble average of the scalar Sk(R) con- 
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structed for the chosen center out of the two moments 
given above: 



S k (R) 



where 
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(23) 



(24) 



(25) 
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where N denotes the number of centers chosen in car- 
rying over the average. The scalar Sk{R) takes on val- 
ues in the range [0, 1] attaining its maximum value of 
unity when the particles are aligned along a straight 
line passing through the center and zero for a uniform 
distribution of particles around the centre. 

Evidently, using the above statistic we cannot in- 
fer about filamentarity given only a density field since 
we must have information about particle positions in 
order to compute S(R). Consequently, we will not be 
able to evaluate this statistic for the adhesion model. 
(It would certainly be worthwhile to construct a statis- 
tic that would work for densities.) For the rest of the 



els on the left in Fig. ||a, which correspond to earlier 
epochs, show predominantly small scale filamentarity 
while those on the right, which correspond to later 
epochs, show relatively large scale filamentarity. On 
comparing the value of the statistic for different spectra 
but at an epoch when the same scale is going nonlinear 
in all of them we see that S is larger for spectra with 
lower value of n. Moreover, at later epochs there are in 
general two preferred scales at which the filaments occur 
predominantly. This is seen very clearly in the n = — 2 
model (Fig. ||b, bottom panels) wherein we see a lot of 
power on very small scales, which sharply drops close 
to zero on intermediate scales but rises again showing 
substantial filamentarity on large scales. In short, the 
statistic generically shows, at late times, two prominent 
peaks (see right panels of Fig. ||a and ||b). This is en- 
tirely consistent with what is expected from Fig. ^|a-d. 
This demonstrates that the statistic S truly character- 
izes filamentarity. (One might want to increase the effec- 
tiveness of this statistic by first using an algorithm such 
as the Minimal Spanning Tree to delineate the main fea- 
tures of the distribution, and then apply the filamentary 



statistic exclusively to these features (Pearson fc Coles 



approximations and N-body we have chosen, at each 1994: Salmi & Coles 1994) 



epoch, a random sample of about 2% of all the parti- 
cle positions as "test" centers in computing S(R). The 
behaviour of S(R) for N-body, FF, LP and TZ simu- 
lations is shown in Fig. || for two different epochs as 
indicated by the value of a quoted within the panels. 
The two epochs chosen are exactly as in Fig. || (see the 
caption of Fig. |J). The scale R essentially characterizes 
the length, in grid units, of the filaments being explored. 

The N-body curves (thick solid lines) in Fig. || show 
a clear transfer of power from smaller to larger fila- 
ments as successively larger scales go nonlinear. This 
behaviour, reflected by the evolution of the filamen- 
tary statistic, is consistent with the visual impression 
rendered by Fig. ||a-d. For instance, the N-body pan- 



At earlier epochs both FF and LP predict the right 
value of the statistic on most scales and for all spec- 
tra. In the case of n = — 2 models there is an excellent 
agreement, at all times, between FF, LP and N-body 
simulations. However, at later epochs none of the ap- 
proximations produce the right behaviour of the statis- 
tic: FF predicts more power on smaller scales and lesser 
power on larger scales while LP predicts lower power 
on almost all scales. The reason why FF produces more 
power on smaller scales can be seen from panels cor- 
responding to FF in Fig. ||a-d: FF simulations show a 
lot of filamentarity on small scales and none on large 
scales. Due to an inherent assumption of this approxi- 
mation there cannot be any shell crossing. This means 
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that particles approach the "valleys" of the potential 
with ever decreasing speeds leading to many whisker 



is quite sensitive to the primordial spectrum of den- 



sity fluctuations (Kauffmann & Melott 1992; Sahni 



shaped objects. Matter that gets collected along the Sathyaprakash & Shandarin 1994). Hence a proper un 



streamlines of the velocity field in this way eventually 
drains out into the minima of the potential wells. As 
a result, large scale filaments never mature in this ap- 
proximation. The problem with FF is just the opposite 
to that with the Zel'dovich approximation: In the latter 
the pancakes do not remain thin and in the former the 
pancakes never thicken. Both these aspects are contrary 
to the findings of N-body simulations. In fact, the fun- 
damental reason why LP and FF cannot produce large 
scale filaments is because they never use power on large 
scales and hence cannot build up matter coherently on 
such scales. 

TZ predicts a much lower value of the statistic on 
all scales and at all epochs. One may think this contrary 
to what the figures |^a-d seem to indicate, especially on 
very large scales. However, we recall that in TZ sim- 
ulations the amount of matter in overdense regions is 
under 30% at most times (cf. Fig. 0) indicating that a 
bulk of the matter is distributed diffusely in undcrdense 
regions. This means that the strong visual signal that 
we see (especially in Fig.||c) is not picked up by this 
statistical test. Consequently, we obtain, in this case, 
a very low value of S. (One might anticipate a better 
agreement between TZ and N-body if the low signal- 
to-noise ratio in TZ were amplified using the Minimal 
Spanning Tree or a fixed overdensity threshold on which 
to test for filamentarity.) 

3.2.5 Evolution of voids 

The next in our list of structural units are voids. There 
is now a consensus on the view that most of the uni- 
verse is filled with voids of different shapes and sizes, 
with most of the matter residing on boundaries sepa- 
rating them. The distribution of volume amongst voids 
of various sizes, sometimes called the "void spectrum" , 



derstanding of the distribution of voids in the Universe 
could in principle lead to a precise estimation of the pri- 
mordial spectrum. Further, Sahni, Sathyaprakash and 
Shandarin 1994, have pointed out that the void spec- 
trum can potentially be used to determine the value 
of the density parameter f2 if the shape and amplitude 
of the initial perturbation spectrum are independently 
known. Thus, the statistics of voids is an important 
measure in characterizing the large scale structure of 
the Universe. In our study of voids we employ two indi- 
cators of void statistics: (i) The total number of voids 
in our simulations (which is a function of epoch), and 
(ii) the void probability function. The former will be 
discussed here and the latter in the next subsection. 

We define individual voids as connected regions of a 
given underdensity. As in the case of clumps the number 
of voids is sensitive to the threshold density chosen for 
identifying them. Again, since the aim here is to com- 
pare the various approximation schemes with N-body, 
we need only choose an appropriate threshold density so 
that voids are (visually) well defined. We found that in 
order to obtain a good picture of voids it is necessary to 
smoothen the density field before applying the density 
threshold in selecting void regions. If the density field 
were not smoothed then there would be too many tiny 
voids which would give rise to a lot of "noise" in the evo- 
lution of the number of voids without at the same time 
making any significant contribution to the total volume 
occupied by voids. In Fig. || we have shown regions in 
our N-body simulations of density p < p c at two differ- 
ent epochs each for the n — 0, k c = 32k / model (top 
panels) and n = 0, k c = 256k t model (bottom panels). 
The density fields were smoothed by removing power on 
modes k > 64k f. We see that voids are well defined for 
the chosen underdensity except that there are still a few 
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voids of very small size. Thus, in our definition of voids 
we do not take into account voids whose diameters are 
smaller than 10 grid units. 

The evolution of the number of voids plotted in 



Fig. 10 exhibits behaviour similar to the evolution of 



clumps seen in Fig. [?]: With reference to the N-body 
curve (thick solid line) we observe that initially there is 
a sharp increase in the number of voids, it acquires a 
maximum sometime after the epoch <r* (except in the 
n = — 2 model where there is no peak), and thereafter 
falls steadily, more or less stabilizing after a while. The 
voids are not well defined at very early epochs (<r < 1) 
but by the epoch <r* they gain their identity. The decline 



in the number of voids in Fig. 10 is a consequence of the 
fact that voids compete for space during expansion, and 
that smaller voids can be encroached upon by larger 
ones. Thus, voids not only expand they can also con- 



tract and ultimately disappear (Bahni, Sathyaprakash 



& Shandarin 1994) 



The adhesion model produces roughly the right evo- 
lution for the number of voids, the largest disagreement 
being for the n = 0, k c — 256k f model. In the latter 
case it predicts a slightly smaller number of voids than 
what is predicted by N-body simulations. FF and LP 
fail to reproduce the correct evolution of the number 
voids in all but the n — 0, k c — 256k f model where 
they agree with N-body very well. In the case of FF 
as matter falls into deeper wells of the potential, the 
cellular structure gets completely phased out with the 
result that at later epochs very few voids are left behind 
in this case. In other words, in FF there is no cellular 
structure to "support" the voids. The same is true in 
the case of LP but because of the fact that the particles 
do cross over caustic regions the cellular structure lasts 
a little longer and the fall off in the number of voids is 
somewhat slower than in the case of FF. Note especially 
that in the n = 2 model both FF and LP predict only 
one void at all epochs. This is consistent with the visual 



pictures corresponding to these approximations (cf. Fig. 
|2|c) which indicate that these simulations will only have 
one void with a sponge like topology. TZ, as expected 
from the pictures in Fig. |^, always predicts fewer but 
larger voids as compared to N-body. This is because 
regions that are populated by many small clumps in 
N-body are filled with a low rather uniform density in 
TZ. It is therefore not a suitable tool to study the void 
spectrum. 

Summarising, the adhesion model is best suited for 
studying the statistics of voids: It not only predicts the 
right evolution for the number of voids (the sharp rise 
and the subsequent gradual fall-off ) , it also predicts the 
right number of voids for most spectra at virtually all 
times. 

3.2.6 Void probability function 

At late times voids are large scale coherent features 
and it is hard to draw a definite demarcation bound- 
ary about a void. There is great danger in using the 
"friends-of-friends" algorithm to identify voids, espe- 
cially at late times. A closer look at Fig. || shows that 
a tiny bridge connecting two neighbouring voids will 
cause the algorithm to declare as one void what visually 
would appear to be two distinct voids. For such epochs 
the void probability function (henceforth referred to as 
VPF) is a better indicator of the sizes of voids than 
void number. We therefore supplement the information 
obtained by studying the evolution of the number of 
voids with that obtained using the VPF. The VPF de- 
scribes the probability that a sphere of size R (a circle 
in 2D) thrown at random is completely devoid of matter 



(White 1979). In practice one can relax this condition a 
little and say that the VPF is the probability of finding 
that a sphere of radius R placed at random within the 
simulation box is an "underdense" region. Here again 
we are faced with a non-objective definition since the 
results will depend upon the chosen value of underden- 
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sity. However, for our purposes of comparison, it makes 
sense to choose any reasonable underdensity that would 
be consistent with what the pictures project. We first 
make a map of the overdense regions as in, say, Fig. |^ 
and for this map we compute the VPF. We have chosen 
at random 20 % of all the grid points in computing the 
VPF. 

Our results are shown in Fig. [ll| a and b (two epochs 
each) for all the models considered by us. The epochs 
chosen are the same as in Fig. ||. We notice that the VPF 
at earlier epochs falls off sharply with scale indicating a 
scarcity of large voids. At later epochs however, the fall 
off is slower, indicating the formation of larger voids as 
the Universe expands. TZ predicts, as expected, a lot 
more voids on all scales and at all epochs as compared 
to N-body. The reason for this is simple: At any epoch 
a all the power in the primordial potential below the 
scale A:nl(c) _1 has been removed in evolving the par- 
ticles up to the epoch a and hence one cannot expect 
the formation of structure on this and lower scales in 
TZ simulations. In the n = models, at earlier epochs 
AM, FF and LP, all agree with N-body. However, at 
later epochs the latter two approximations predict fewer 
voids of all sizes. The reason for this is that matter does 
not participate in hierarchical clustering in FF and LP 
and thus matter does not get emptied out as happens in 
the case of N-body. AM comes closest to the predictions 
of N-body for this spectrum. 

While this is the story in the case of spectra with 
equal power on all scales the results are totally different 
when the power is not equal on all scales. In the case of 
the n = 2 model none of the approximations agree with 
N-body with AM and TZ being closest in accuracy. For 
n = — 2 on the other hand FF and LP agree remarkably 
well with the predictions of N-body both at earlier and 
later epochs and on almost all scales. 



4 CONCLUSIONS 

In this work we have studied gravitational instabil- 
ity in the strongly non-linear regime using a number 
of distinct and sometimes orthogonal statistical indi- 
cators such as: Correlation coefficient of particle posi- 
tions/densities; statistics of overdense and underdense 
regions (clumps and voids); filamentary statistics; and 
the Void probability function. Using these statistics we 
assess the accuracy of different approximations to grav- 
itational instability such as: the Truncated Zel'dovich 
approximation (TZ), the adhesion model (AM), the 
frozen flow (FF) and the linear potential (LP) approx- 
imation. We compare these approximations with N- 
body simulations for a variety of spectra and at dif- 
ferent cosmological epochs. We find that as the scale of 
non-linearity grows, so do the characteristic features of 
structure such as filamentarity and the sizes of voids. 
Our study shows that as long as the non-linear length 
scale fcj^L remains smaller than - the typical dis- 
tance between peaks of the gravitational potential - 
all approximations give results in reasonable agreement 
with those of N-body simulations. During the epoch 
i?* < fcj^L < R v (R v being the correlation length of the 
potential) the adhesion model (and occasionally TZ) 
gives results closest to N-body. The reason for this is 
the following: Particles in FF follow the streamlines 
of the initial velocity field, converging in the minima 
of the gravitational potential. As a result, filaments in 
this approximation, show a tendency to thin out and 
ultimately disappear. In the case of LP, particles begin 
to oscillate about the minima of the potential at late 
times, freezing the possibility of any long range dynam- 
ics. Consequently, no real evolution of particle positions 
can take place beyond the epoch fe^ L ~ R* in either FF 
or LP which break down when > i?,* . 

Unlike FF and LP both AM and TZ use the pres- 
ence of long range modes in the gravitational potential 
to move particles at late times. Consequently, as long 



as the potential has sufficient long range power to af- 
fect bulk motion over large scales (described by Rp), 
TZ and AM give results closely matching with those of 
N-body. Compared with AM, TZ suffers from one ma- 
jor drawback: although it does manage to collect the 
right amount of matter into large clumps it completely 
overlooks the presence of small clumps. The reason for 
this is clear, in order to prevent shell crossing all modes 
which have gone non-linear by a given epoch are sur- 
gically removed from the initial gravitational potential 
because of which no small scale clustering is present 
in TZ. The adhesion model does not suffer from this 
drawback and can accurately predict the multiplicity 
function of clumps and the void spectrum as long as 
— R<p- However, after fc^L — ^-v TZ makes much 
more accurate predictions about the location of mass. 
In this regime, AM continues to make better statistical 
predictions, but its dynamical accuracy is reduced. This 
can be understood as inaccurate influence of the short 
modes on the position of structure in AM at late times 
especially since by construction, the adhesion technique 
when applied to truncated initial spectrum closely re- 
produces TZ. 

A comparison of i?„ and R v for different spectra 
(Fig. [l]), shows that i?* ~ R v for very steep or very 
shallow spectra with n > 2 or n < - 2. For interme- 
diate values — 2 < n < 2 R v ^> i?* indicating that 
for such spectra the adhesion model (and occasionally 
TZ) will be more accurate than FF or LP. These val- 
ues of the two dimensional spectral index correspond 
in three dimensions to the range —3 < n < 1 which is 
precisely the range of interest in most cosmological sce- 
nario's such as CDM. (For the standard CDM model, 
i?* < 1 Mpc. R v ~ 50 Mpc.) We therefore feel that the 
adhesion model and the truncated Zel'dovich approxi- 
mation (depending on which aspect of the description 
of structure one wishes to emphasize) are more realis- 



tic approximations to apply to the study of large scale 
structure than either FF or LP. 
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Figure 1. The two natural scales of the potential _R* and R v (left 
panel) and the "temperature" 7 (right panel) are shown plotted 
against the spectral index n for featureless power-law spectra. 
The values of ij* and R v are quoted for a box of size 512 X 512. 
Note that the difference between the two scales decreases as |n| 
increases, being largest when n ~ (i?» and R v characterize the 
mean distance between peaks and the correlation length of the 
potential, respectively.) 



Figure 2. Comparison of N-body simulations with the simula- 
tions of various approximation schemes. Only a fourth of all the 
particles are shown. For the sake of clarity, we have superim- 
posed the results of adhesion model over the particle positions 
of N-body simulations. From top to bottom the pictures corre- 
spond to N-body, adhesion model, frozen flow, linear potential, 
and truncated Zel'dovich approximation, respectively. The left 
and the right panels are obtained at epochs cti and <T2, respec- 
tively. The plots are shown for (a) n = 0, k c = 32k f, 01 = 4, 
0-2 = 16, (b) n = 0, k c = 256fc/, a 1 = 32, <t 2 = 128, (c) n = 2, 
k c = 256k f , cti = 256, <7 2 = 1024,and (d) n = -2, k c = 256fc/, 
cti = 2.00, (7 2 = 2.83. 



Figure 3. Contour plots of the initial potentials corresponding to 
the different initial potentials used in our simulations. The solid 
lines correspond to ip > and the dashed lines to <p < 0. 



Figure 4. Time evolution of the vector correlation coefficient of 
particle positions corresponding to (i) frozen flow (dashed line), 
(ii) linear potential (dotted line), and (iii) truncated Zel'dovich 
(dashed-dotted line). Top left panel corresponds to n = 0, k c = 
32k f power law model and top right panel to n = 0, k c = 256fc f 
power law model. Bottom left panel is for n = 2, k c = 32k f and 
the one on the bottom right is for n = —2, fc c = 256k f. The 
dotted and the solid vertical lines correspond to the epochs when 
the scales going nonlinear are R„ and R v , respectively. 



Figure 5. Evolution of the density correlation coefficient cor- 
responding to (i) adhesion model (solid line), (ii) frozen flow 
(dashed line), (iii) linear potential (dotted line), and (iv) Zel- 
dovich (dashed-dotted line). In Fig. 5a the panels arc arranged as 
in Fig. 4 and before the correlation coefficient is computed the 
n = 0, fc c = 32fc^ density field in each simulation is smoothed 
at fcc = 16fcy and the rest of the density fields are smoothed at 
kc = 64fc^. In Fig. 5b we have shown the evolution of the density 
correlation coefficient, for the n = 0, fc c = 256fc f model, obtained 
by smoothing the density fields, at a given epcoh, at a scale pro- 
portional to the nonlinear scale at that epoch : fc^ 1 = 0.5 X fcj^- 
This demonstrates that while the TZ reproduces the large scale 
features very accurately, the other approximations, with the ex- 
ception of the AM, even after such a smoothing, do not show good 
agreement with N-body results. 



Figure 6. Plot showing regions of density greater than the 
threshold density p c in our N-body simulation for n = mod- 
els. Notice that for the chosen threshold density clumps are well 
defined features. Top panels correspond to fc c = 32k f and bottom 
panels to k c = 256fc f. 
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Figure 7. Evolution of the number of clumps (top panels) and 
the total mass in clumps (bottom panels) corresponding to (i) 
N-body simulations (thick solid line), (ii) adhesion model (solid 
line), (iii) frozen flow (dashed line), (iv) linear potential (dotted 
line), and (v) truncated Zcl'dovich (dashed-dotted line). The evo- 
lution is shown for: (a) n = spectrum with left and right panels 
corresponding, respectively, to, k c = 32k f and k c = 256fc f, and 
(b) for n = 2, fc c = 256k f spectrum (left panels) and n = — 2, 
k c = 256k f spectrum (right panels). 



Figure 8. Filamentary statistic S(R) is shown plotted as a func- 
tion of the scale length R for two different epochs (as indicated 
by the value of a within the panels) corresponding to (i) N-body 
(thick solid line), (ii) frozen flow (dashed line), (iii) linear poten- 
tial (dotted line), and (iv) truncated Zel'dovich (dashed-dotted 
line). Results are shown for (a) n = with k c = 32k f (top pan- 
els), fc c = 256k f (bottom panels) and (b) n = 2, k c = 256k f 
spectrum (top panels) and n = —2, k c = 256k f spectrum (bot- 
tom panels). 



Figure 9. Plot showing regions in our N-body simulations of 
density less than a threshold density p c for n = models at two 
epochs each for cutoff at k c = 32k j (top panels) and k c = 256k j 
(bottom panels). 



Figure 10. Evolution of the number of voids corresponding to (i) 
N-body simulations (thick solid line), (ii) adhesion model (solid 
line), (iii) frozen flow (dashed line), (iv) linear potential (dot- 
ted line), and (v) truncated Zel'dovich (dashed-dotted line). The 
panels are as in Fig. 4. 



Figure 11. Void probability function V(R) is shown plotted as 
a function of the scale R for two different epochs (as indicated 
by the value of a within the panels) corresponding to: (i) N-body 
simulations (thick solid line), (ii) adhesion model (solid line), (iii) 
frozen flow (dashed line), (iv) linear potential (dotted line), and 
(v) truncated Zel'dovich (dashed-dotted line). Results are shown 
for (a) n = with k c = 32fcy (top panel), k c = 256k f (bottom 
panel) and (b) n = 2, k c = 256k j spectrum (top panel) and 
n = —2, k c = 256kf spectrum (bottom panel). The overdensity 
field needed in finding VPF is determined by taking a threshold 
density of p c = 5po in all but n = 0, k c = 256k f model for which 
the threshold is taken to be p c = 2po- 



